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We present numerical results from three-dimensional evolutions of scalar perturbations of Kerr 
black holes. Our simulations make use of a high-order accurate multi-block code which naturally 
allows for fixed adaptivity and smooth inner (excision) and outer boundaries. We focus on the 
quasinormal ringing phase, presenting a systematic method for extraction of the quasinormal mode 
frequencies and amplitudes and comparing our results against perturbation theory. 

The detection of a single mode in a ringdown waveform allows for a measurement of the mass 
and spin of a black hole; a multimode detection would allow a test of the Kerr nature of the source. 
Since the possibility of a multimode detection depends on the relative mode amplitude, we study this 
topic in some detail. The amplitude of each mode depends exponentially on the starting time of the 
quasinormal regime, which is not defined unambiguously. We show that this time-shift problem can 
be circumvented by looking at appropriately chosen relative mode amplitudes. From our simulations 
we extract the quasinormal frequencies and the relative and absolute amplitudes of corotating and 
counterrotating modes (including overtones in the corotating case). We study the dependence of 
these amplitudes on the shape of the initial perturbation, the angular dependence of the mode and 
the black hole spin, comparing against results from perturbation theory in the so-called asymptotic 
approximation. We also compare the quasinormal frequencies from our numerical simulations with 
predictions from perturbation theory, finding excellent agreement. For rapidly rotating black holes 
(of spin j = 0.98) we can extract the quasinormal frequencies of not only the fundamental mode, 
but also of the first two overtones. Finally we study under what conditions the relative amplitude 
between given pairs of modes gets maximally excited and present a quantitative analysis of rotational 
mode-mode coupling. The main conclusions and techniques of our analysis are quite general and, as 
such, should be of interest in the study of ringdown gravitational waves produced by astrophysical 
gravitational wave sources. 

PACS numbers: 04.30.Db, 04.70.-s, 04.80.Nn, 04.25. Dm 



One of the most useful methods to explore the response of black holes to external perturbations is based on wave 
scattering 0. Early studies identified three main stages in the dynamics of a wave propagating on a black hole 
background, as observed at a fixed spatial point. In a first, transient phase the observed wave depends on the 
structure of the initial pulse. Vishveshwara and Press discovered that this initial "burst" is invariably followed by a 
second phase characterized by exponentially decaying oscillations: this phase is usually referred to as "quasinormal 
ringing" |^ . In the third and last stage of the evolution, waves slowly die off as a power law tail Q . 

Astrophysical black holes should be well described by the Kerr solution, since charge is unlikely to play a major 
role in astrophysical scenarios (sec e.g. for a discussion). As a consequence of the "no hair theorem", if general 
relativity is the correct theory of gravity, the quasinormal mode (QNM) frequencies of a Kerr black hole depend only 
on its mass and angular momentum. Earth-based and space-based gravitational wave detectors have the potential to 
measure the frequency and damping time of a QNM. From these two observables we can infer the black hole's mass 
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and angular momentum la,l3- Fo'" the space-based Laser Interferometer Space Antenna (LISA) , and possibly also 
for second-generation ground-based detectors, the signal-to-noise ratio can be large enough that we will be able to 
identify two or more QNM frequencies in the signal A multi-mode detection would provide a striking, direct test 
of the Kerr nature of the source (i.e., of the no-hair theorem). The basic idea is quite simple. Roughly speaking, the 
first mode in the pair is used to determine the black hole's mass and angular momentum, and the other mode(s) to 
verify that the QNM spectrum is indeed consistent with a general relativistic Kerr black hole 0. 

The combined observation of supermassivc black hole binary inspiral and ringdown with LISA can provide even 
more information p^ . Parameter estimation during the inspiral phase can be very accurate, depending on the black 
holes' masses, spins and distance [Tll |. Combining information from the inspiral and ringdown phases we can estimate 
the energy radiated in the merger, and possibly improve parameter estimation from both phases (see e.g. |l2l | for a 
preliminary study of this effect in the context of earth-based detectors) . 

In the last thirty years the development of gravitational wave astronomy motivated a detailed investigation of the 
QNM frequency spectrum hj. In comparison, the problem of the relative excitation of QNMs received very 

little attention (see e.g. and references therein). Ideally, the relative QNM excitation should be determined 

by general relativistic simulations of binary black hole mergers. Despite recent progress, this information is not yet 
available [T^ . Given the recent progress of numerical relativity, by the time LISA flics we could have a good knowledge 
of the multipolar distribution of the energy and angular momentum radiated in a black hole merger under generic 
conditions. Knowing in advance which modes should be excited in a realistic merger will not only be useful to probe 
the Kerr nature of the source, but also to reduce the number of templates needed to perform matched filtering on 
ringdown waveforms. 

In this paper we present a quantitative investigation of QNM excitation studying a simple model problem: the 
scattering of scalar waves on a Kerr background. We use our new infrastructure for multi-block simulations |l8l | for 
these studies. The infrastructure is based on the techniques described and applied in the context of numerical relativity 
in 10 and further extended in . Our infrastructure uses Carpet 0, |23| , a driver for the Cactus computational 
toolkit US 13' originally designed to provide fixed and adaptive mesh refinement. The capabilities of Carpet have 
been recently extended [l^ to include the type of boundary conditions that are needed for multi-block (also called 
multi-patch) simulations in Cactus. 

Multi-block techniques yield increased efficiency and accuracy in our studies, for two main reasons. The first is 
that we can set up smooth excision and outer boundaries, and we can therefore apply boundary conditions in a clean 
and well understood way. Furthermore, a suitable multi-block grid structure provides a natural and flexible way of 
implementing mesh refinement. We can keep a fixed angular resolution throughout the entire domain, avoiding the 
unnecessary high resolution at large distances from the central object that one would have using a cartesian grid. 
The resources thus saved can be used, for example, to set up a rather large number of grid points in the radial 
direction, so that outer boundaries are located at large radii and the noise produced at the boundaries does not affect 
the results. Even though we are working in three dimensions, our results are more accurate than previous studies 
using two-dimensional codes . Krivan et al. studied the late time dynamics and the rotational coupling 

of masslcss scalar fields in a Kerr background, but not their quasinormal ringing. Later they extended the analysis 
to gravitational perturbations, considering both the late time tail and the quasinormal ringing phase [26l |. For large 
rotation the damping times of corotating fundamental modes in |26l| are accurate within ^ 3% when compared to 
results from perturbation theory; our accuracy (~ 0.3%) is roughly an order of magnitude better. In fact, we can 
extract the frequencies of some overtones with an error of the order of a few percent or less. 

Given the high accuracy of our multi-block infrastructure, a careful extraction of the QNM content of the waveforms 
becomes necessary. We discuss in detail the so-called time-shift problem (exponential dependence of the quasinormal 
amplitudes on the time at which the quasinormal ringing regime starts), how it affects the determination of both 
absolute and relative QNM amplitudes, and how to choose pairs of modes so as to decrease the uncertainty on relative 
amplitudes. We also introduce a general criterion (based on minimizing a suitably defined residual) to determine the 
optimal fitting window to extract QNM frequencies and amplitudes. Using these tools we study the absolute and 
relative amplitudes of corotating and counterrotating modes for Gaussian initial data located in the far zone. We 
study the dependence of these arnplitudes on the radial shape of the initial data, finding excellent agreement with 
results from perturbation theory [l6j. We also discuss the problem of extracting overtones for modes with a given 
angular dependence, finding that the first overtones of corotating modes (e.g. modes with / = m = 2) contribute 
significantly to the waveform for rapidly rotating black holes. 

The plan of the paper is as follows. In Sec. we briefly describe our multi-block code. After introducing the 
background metric, we discuss the numerical implementation of the scalar wave equation and our time evolution 
techniques. In Sec. IIIII we extract QNM frequencies and amplitudes from our evolutions, comparing with analytical 
predictions from perturbation theory. To start with, we point out some conceptual limitations in the extraction of 
QNM amplitudes due to the so-called time-shift problem. Then we introduce a rather general method to determine 
the best fitting interval to extract QNM waveforms. We first check the accuracy of this method (and of our numerical 
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code) by reproducing the QNM frequencies predicted by standard perturbation theory. Scalar QNM frequencies for 
Kerr black holes have been computed in Q, and they have never been systematically confirmed by numerical time 
evolutions^. Next we give a quantitative estimate of rotational mode mixing as a function of the black hole's spin and 
discuss the initial data dependence of the amplitudes of corotating and counterrotating modes. In this way we assess 
the validity of the amplitudes predicted by perturbation theory in the so-called asymptotic approximation (where 
both the observer and the initial data are located far away from the black hole). Finally we discuss the extraction of 
overtones from our waveforms. 



II. NUMERICAL EVOLUTION OF SCALAR PERTURBATIONS OF KERR BLACK HOLES 

A. Grid structure 

We perform our evolutions describing scalar perturbations of a Kerr spacetime through excision of the singularity. 
With our multi-block approach we can have smooth (in particular, spherical) inner (excision) and outer boundaries. 
As in 0, we use a six-block setup with a global topology of x referred to as cubed sphere coordinates (Fig.^). 
This topology and the corresponding coordinates on each block are well adapted for modeling a single central object 
together with outgoing radiation that is generated at or close to that object. 

The six blocks are arranged like the six faces of a cube, i.e., block covers the neighborhood of positive x, block 1 
positive y, block 2 negative x, block 3 negative y, block 4 positive z, and block 5 negative z. On each of those blocks 
a local coordinate system (a,6, c) is defined, with —1 < (a,5, c) < +1, and equal grid spacing in the local system. 
The coordinate c runs along the radial direction, and a, b span the angular ones. See |l9l | for the explicit definition of 
these coordinates. 




Figure 1: Illustration of the six-block grid structure and the cubed sphere coordinates that are used for the simulations in this 
paper. The left panel shows the distribution of grid points on a sphere of constant radius. The central panel shows a snapshot 
from a scalar wave evolution on an equatorial cut. The plot refers to an ^ = m = 2 mode on the background of a Kerr black 
hole with spin j = 0.9 at t = 92. 2M. Also shown are the locations of the inter-block boundaries. The right panel magnifies the 
central region of the domain in the equatorial plane, showing the grid structure around the spherical excision boundary. The 
four dark lines mark the interfaces between blocks. 



B. Background metric 

We consider a stationary, rotating black hole background. The Kerr metric can be written in Kerr-Schild form as 

ds^ = Vt.v + 2Hl^l^dx>'dx'' (1) 



^ See however l27l , where the fundamental scalar mode with / = was observed to dominate the emission of scalar radiation by perturbed 
Kerr black holes in the Brans-Dicke theory of gravity. 
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with rj^i, the Minkowski metric, and 



(2) 

+ {z/r) 



= l{p'-a') + ^^ip'-a^r + a^z\ (3) 

2 2 I 2 I 2 

p = s + y + z . (4) 

Here A/ is the mass and a — jAI ~ J/M is the angular momentum per unit mass of the black hole (j is the 
dimensionless spin parameter, < j < 1). In Cartesian coordinates, the null vector is given by 

. , , , rx + ay , ry — ax , z , 

l^dxt" ^dt+ — l-dx + -dy + -dz . (5) 

This form of the Kerr-Schild metric has become of common use in numerical relativity. However, in these coordinates 
the shape of the Cauchy and event horizons become more and more ellipsoidal with increasing spin^. For j > 0.96 
it is not possible to fit a spherical excision boundary between these horizons any more. This is illustrated in Fig. |21 
Although we could in principle choose a different shape for the excision boundary within our code, we instead use 
coordinates in which both horizons are always spherical, and therefore an excision sphere can always fit between 
them. This version of the Kerr Schild coordinates is related to the "standard" one defined above by the following 
transformation: 

X = X ^ , (6a) 

r 

y = y H , (6b) 

r 

z = z . (6c) 



Event horizon - 
Cauchy horizon - 
Excision boundary 




Figure 2: Event and Cauchy horizons for a Kerr black hole with spin j > 0.96 in "standard" Kerr-Schild coordinates (as defined 
in the text), here shown in the x-z plane. The horizons have an ellipsoidal shape; it is therefore not possible to fit a spherical 
excision region (dotted line) between the two horizons. 



C. Evolution system 



We write the time evolution equations for scalar perturbations in a symmetric hyperbolic and (in the case of 
a stationary background) conservative form. This guarantees stability and energy conservation for the continuum 
equations. We use differencing operators that satisfy the summation by parts (SEP) property, and this also guarantees 
stability and energy conservation in the semi-discrete case (see e.g. |28l | for more details). On a time independent 
background the evolution equations take the form 

$ = H, (7) 

ri = /3*a,H+-^9, ( —/3'U + aVhH'^dj] , (8) 
Vh \ a J 

d, = a,n, (9) 



^ We thank Harald PfcifTcr for pointing this out to us. 
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where $ denotes the scalar field, 11 its time derivative, and di = 9^$ the spatial gradient of the field. The quantity 
hij is the three metric, h its determinant, h^^ the inverse three metric, a the lapse, and /3' the shift vector. H^^ = 
fiV _ p'^pi ja^ is the spatial part of the inverse four-metric. 

The background geometry for all simulations presented here is that of a Kerr black hole. It would be possible to 
exploit the axisymmetry of the background spacetime by performing a multipole decomposition of the scalar field, 
and then solving for each azimuthal number m as an axisymmetric, two-dimensional problem. We choose not to do 
so here but instead solve the full three-dimensional equations. This has the advantage that wc can later use the 
same implementation for generic, non-axisymmetric spacetimcs. Using a fully three-dimensional code also serves to 
test our numerical multi-block and excision techniques in a scenario that is non-trivial, but at the same time not as 
complicated as solving the full Einstein equations. 



The QNM excitation depends on the angular structure of the scalar field that is used as a perturbation. To excite 
certain modes in a controlled way, we choose initial data of the form 



Unless otherwise stated, throughout this paper we use rg = 20M and a = M. Y(rn{S, <t>) denotes the ordinary spherical 
harmonics. Since the Kerr background is not spherically symmetric, we should really expand the perturbation in terms 
of spin-weighted spheroidal harmonics sSim{j^) of spin weight s = 0. Using spherical harmonics weakly excites other 
modes through rotational mode mixing; this point will be discussed in more detail below, in Sec. IIII DI 

The changes in the characteristic length scale in the radial direction are usually small over time. To accurately 
resolve the propagating waves all the way to the outer boundary we use a constant resolution in the radial direction 
of our cubed sphere coordinates. As mentioned, the coordinates are set up so that the spherical inner (excision) 
boundary is placed between the event and Cauchy horizons, and no boundary conditions need to be applied there. 
For global stability wc choose maximally dissipativc boundary conditions at the outer boundary, and we apply them 
through penalty terms. 



We use spatial finite differencing operators that satisfy summation by parts; they are eighth order accurate in the 
interior and fourth order accurate at and close to the boundaries. With those operators we expect a global accuracy 
of order five (see [23| for more details on the operators that we use). We use a fourth order accurate Runge Kutta 
time integrator. This does not spoil the expected global fifth order spatial convergence, since we use a small enough 
time step so that the truncation errors generated by the time integration are smaller than those that originate from 
the spatial finite differencing (see |23| for details on the code's convergence). 

In multi-block simulations one does not necessarily have a uniform or isotropic grid spacing in a global coordinate 
system. Since in all our simulations the global grid spacing in the radial direction is smaller than in the angular 
directions, we use the radial direction for our time step criterion Ai = AAr, where A — usually referred to as the 
Courant factor — is chosen to be A = 0.25. 

Taking into account our initial data [cf. Eq. (|10|l ] and the typical position of the observer, and given that we are 
interested in the ringdown phase, evolution times of about t = 150A/ with outer boundaries at about 200Af are a 
reasonable choice. Unless otherwise stated, for the simulations that wc show below wc use ten points per length unit 
M in the radial direction and an angular resolution of 21 x 21 grid points per block, which gives us approximately 80 
grid points along the circumference of any sphere with constant radius. As described below, we have found that with 
this resolution we can get good agreement with the Kerr quasinormal frequencies predicted by perturbation theory. 

Figure O shows a typical waveform that we get when extracting the real part of the £ = 2, 77i = 2 mode from 
our simulations. The initial data are set up according to Eq. H10|l . with the specific choice A — B = 1, a = M 
and ro = 20il/. The background Kerr black hole has a spin j ~ 0.9. The strongest modes in this waveform are 
(.£ = 2,TO = 2,n = 0) and (£ = 2,m = —2,n = 0), where we use n to label overtones, n = being the fundamental 
mode. We show a fit for those two modes together with the numerical data. The third strongest component in the 



D. Initial and boundary conditions 




(10b) 



(10a) 



(10c) 



E. Specifications for the simulations 



6 



data is the = 2, m = 2, n = 1) mode. Since this mode is decaying much faster than the fundamental mode, it only 
plays a role at early times. That is the reason why our fit, done for only the two fundamental modes, is drifting away 
from the numerical data at times below 50A/ (we will explicitly analyze overtones in Sec. IIII F|) . 




Figure 3: The left panel shows the I = m = 2 component of the waveform extracted at radius r — 5M on a Kerr black hole 
background with a spin of j — 0.9. The waveform is a superposition of the corotating and the counterrotating mode, and the 
beating of two different frequencies is clearly visibie. The right panef shows the waveform for t > AOM as wefl as a QNM fit 
with the fundamental £ = |m| = 2 modes. The interval used for the fit is [74.5Af, 150Af]. The inlay shows the absoiute value 
of the difference between the fit and the data. At times between the excitation of the QNM (t ~ 25M) and about 70M the 
differences are mainly due to the presence of the {I — 2,m = 2,n = 1) mode, the exponentially decaying mode that can be seen 
in the inlay (a fit of this mode yields quasinormal frequencies in agreement with perturbation theory). At times t < 25M the 
difference is due to the initial burst. 



III. QUASINORMAL MODE FREQUENCIES AND EXCITATION AMPLITUDES FROM NUMERICAL 

SIMULATIONS 



A. Overview 



The time evolution of perturbations of a Kerr black hole can be split into three stages. After a first burst of radiation 
depending on the source of the excitation, the perturbation field $ undergoes exponentially damped oscillations 
(ringdown phase). Finally, in the tail phase (caused by backscattering of radiation off the background gravitational 
potential) the field follows a power-law decay. In this paper we focus on the ringdown stage. We extract the different 
multipole components of the numerical solution by integrating the scalar field against different spherical harmonics 
over surfaces of constant observer radius r: 



,t) = I y,:„(0,0)$df^, (11) 



where a star denotes complex conjugation. We usually consider multipole components up to £ = 4 and all values of 
w (|w| < £). By adding up the contributions of all multipoles one should recover the full scalar field: 



„ oo £ 
•J r »_o a 



As already mentioned, the rotation of the black hole and numerical errors can excite multipole components which 
are not present in the initial data. The above property can be used to check for the existence of overtones or modes 
with € > 4 that are not explicitly extracted but might be present in the solution (e.g. due to rotational mode mixing, 
numerical errors, or both). Multipoles with m 7^ require some care. The spherical harmonics Yfm(^, fp) ^-rc given by 
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where Pp{9) is a real function (an associated Legcndre polynomial). Therefore the initial data of a pure multipole 
with m ^ will be complex. Given that the evolution equations are linear, we can evolve the real and imaginary parts 
of $ separately, and obtain the complex solutions for positive and negative m by linear combinations of the form 

<^i„i = 7^($£™)+iX($^,„), (14a) 

= ni^ira) - ^i^ern) ■ (14b) 

This point is important for the extraction of the relative amplitude of corotating and countcrrotating modes. In fact, 
as stressed (for example) in Ref. |^ , QNMs of Kerr black holes always come "in pairs" . In the Kerr case, for a 
given multipole (£, m) we have to solve an eigenvalue problem to determine both the quasinormal frequencies uJimn 
and the angular separation constant Aimn (not to be confused with the mode amplitude Aimn introduced below), 
used to separate the angular and radial dependence of the Teukolsky equation and write it as two ordinary differential 
equations. For each {£, to 7^ 0) and j 7^ the eigenvalue problem admits two sets of solutions. In addition to {£, to), we 
label the modes of each set by the overtone index n, denoting the frequencies by w^^„ (« = 1,2). For given {£, m, n), 
the solutions corresponding to the two different sets have different values of i^imn (and also of Aimn)- 

(1) / (2) 
^irnn r ^ Irnn ' 

Both the real and imaginary parts are different. In fact, the real part of one of the frequencies is positive and the 
other one is negative: 

If we consider instead the frequencies corresponding to the pair (^, — jti), they are related to those of (£, to) by a simple 
symmetry property: 

-7^(a;«J=7^(a;ii)„J, I(-IL) = X(c.liLj , (41)*=^^--' = 1> 2 ; z ,) • (15) 

In this sense, any solution with positive to is nothing but the "mirror image" of a solution with opposite real part 
and opposite to (see Fig. 6 of for an illustration of this). For to = (or for any value of to in the Schwarzschild 
case) the two "mirror solutions" are degenerate in modulus of the frequency and damping time. However, in general, 
a multipolar component with a given (£, to) will always contain a superposition of at least two different damped 
exponentials. Because of this, it is enough to consider only one frequency for each mode [{£, to) or {£, —to)], since the 
other two frequencies are obtained through this symmetry property; we follow the standard convention of considering, 
for each mode, the frequency with positive real part. Below we will discuss in detail the excitation of these modes, 
extending previous work by Krivan et al. |26l |. 

When the perturbation field is in the quasinormal ringing regime, it can be expanded as a QNM sum of the form 

3>,™(r, <) « 7^ I Aimne'^'-e-'^'-^'-'"^ I , (16) 
U=o J 

where Aimn is the amplitude of the n-th overtone with angular structure given by the pair {£,m), cimn its phase, 
^emn its complex quasinormal frequency and to (which to a first approximation we assume to be the same for all 
modes) marks the time at which the quasinormal regime starts. 

The extraction of gravitational waves from numerical simulations of the full Einstein equations requires the observer 
to be located far away (in the wave zone). For the extraction of QNM frequencies, on the other hand, it is not 
problematic to place the observer close to the black hole, since an observer at any point in the space time is in general 
expected to measure the same frequencies. In fact, a small r is better suited for extracting quasinormal frequencies 
from our simulations simply because outer boundary effects pollute our waveform later, and the ringing regime can 
be observed for a longer time. The availability of a longer ringdown waveform improves the accuracy of the fitting 
procedure that we apply to extract the frequencies. 

The effect of the observer's location on the result is illustrated in Table where we list the frequencies of {£ = 
2, TO = ±2) fundamental modes for a Kerr black hole with spin ) = 0.9 as measured by observers at radii r = 5M , 
20M and 40M. We picked to ^ r + tq in Eq. lfTH|l and A = 0, B = 1 in Eq. ((113) • The results presented in this Table 
are discussed in more detail below (Sec. IIII C|l . Here we simply remark that quasinormal frequencies measured at 
different radii are very close to the analytical predictions, supporting the statement that the observer does not need to 
be far away from the black hole to extract the correct ringdown frequencies. Indeed, for these particular simulations 
the relative error increases with r: the main reason, as explained, is that observers located at large radii see boundary 
effects earlier, so they can only measure a shorter ringdown waveform with respect to observers closer to the black 
hole. 
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Table I: Quasinormal frequencies computed by Leaver's continued fraction method (here labeled "perturb.") and by our time 
domain simulations, with the associated relative differences. We use 21 x 21 points in the angular direction on each block and a 
resolution of M/10 in the radial direction. For j — 0.9 we compare the frequencies as seen by observers located at different radii 
r. Observers at larger radii measure frequencies with larger errors, since boundary effects start to contaminate the waveform 
earlier. 



r 


3 


I, 


m 


^perturb. 


^numerical 


rel. difference (Re,Im) 


5M 


0.0 


2, 





0.48364 


- 0.09676i 


0.48364 - 


0.09676i 


< 10"^ 








0.5 


2, 





0.49196 


- 0.09463i 


0.49190 - 


0.09469i 


4.27 X 10" 


^,6.34 X 10" 


4 




0.5 


2, 


-2 


0.42275 


- 0.09562i 


0.42281 - 


0.09569i 


1.42 X 10" 


^,7.32 X 10" 


4 




0.5 


2, 


2 


0.58599 


- 0.09349i 


0.58589 - 


0.09339i 


1.71 X 10" 


Sl.07 X 10" 


3 




0.9 


2, 





0.51478 


- 0.08641i 


0.51471 - 


0.08646i 


1.36 X 10" 


^,5.79 X 10" 


4 




0.9 


2, 


-2 


0.38780 


- 0.09379i 


0.38781 - 


0.09339i 


2.58 X 10" 


^4.26 X 10" 


3 




0.9 


2, 


2 


0.78164 


- 0.06929i 


0.78144 - 


0.06955i 


2.56 X 10" 


S3.75 X 10" 


3 




0.98 


2, 


2 


0.89802 


- 0.04090i 


0.90940 - 


0.04018i 


1.27 X 10" 


^1.76 X 10" 


2 




0.98 


2, 


-2 


0.38177 


- 0.09338i 


0.38234 - 


0.09743i 


1.49 X 10" 


^4.34 X 10" 


2 


20M 


0.9 


2, 


-2 


0.38780 


- 0.09379i 


0.38694 - 


0.09471i 


2.22 X 10" 


^ 9.81 X 10" 


-3 




0.9 


2, 


2 


0.78164 


- 0.06929i 


0.78244 - 


0.06670i 


1.02 X 10" 


^ 3.74 X 10- 


-2 


40M 


0.9 


2, 


-2 


0.38780 


- 0.09379i 


0.38406 - 


0.09958i 


9.64 X 10" 


^ 6.17 X 10- 


-2 




0.9 


2, 


2 


0.78164 


- 0.06929i 


0.78292 - 


0.06618i 


1.64 X 10" 


^ 4.49 X 10" 


-2 



B. The time shift problem 



Here we discuss the so-called time-shift problem, how it affects the extraction of quasinormal frequencies and 
amplitudes from numerical simulations, and a possible way to address it. Even though in this paper we consider 
scalar perturbations, the discussions of this and other sections apply also to other types of black hole perturbations. 

The standard approach is to choose to in Eq. Hlt)|) using some approximate calculation based, for example, on the 
location of the initial data and the time it would take for initial data to be scattered by the black hole potential and 
reach the observer, usually assuming that perturbations propagate with coordinate speed one (as they would in flat 
spacctime). Criteria like this are well motivated and provide a good guess, but there is still an uncertainty in to- For 
example, the coordinate speed of the perturbation in a curved background in general will not be one. One might 
expect that such a small uncertainty would not influence the extraction of physically relevant quantities. However, 
as we discuss below, this is not the case: there are quantities of interest to gravitational wave detection which have a 
strong dependence on t^. Following the existing literature, we will call this the time-shift problem. 

Suppose the starting time to is subject to an uncertainty (5o. Under a change 



to ^ io + <5o , 

the amplitude and phase of each mode change according to 



(17) 



(18a) 
(18b) 



That is, an uncertainty in to induces a linear uncertainty in the phase, and an exponential uncertainty in the amplitude. 
Fortunately other quantities are largely independent of this uncertainty: for example, the QNM frequencies LOimn are 
unaffected by So- 

How large can we allow this exponential amplification of errors to be? Let us require the amplitude uncertainty 
induced by the starting-time uncertainty Sq to be less than some small number e, that is 



A' 



A 



-5oI(wf,„„) 



1 



< e. 



For small e this implies 



l^ol 



< 



I(Ma;£„„) 



M . 



(19) 
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For the ^ = 2 fundamental scalar mode in the Schwarzschild background (which is spherically symmetric, so that 
the choice of m becomes irrelevant) \T {Muj2oo) I = 0.09676 ~ 10~^. In other words, if we want to determine the 
amplitude of this mode within 1% (e = 10^^) we need to know Iq with an uncertainty 6q < O.IM . Constraints on Sq 
are even tighter for overtones, since they decay faster and the exponential propagation of errors is more dramatic. 

In practice, what is most interesting is the relative amplitude between different modes. Under a change of the form 
()17|l this relative amplitude changes according to 
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£nin 
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ijnn 



Af 
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(20) 



Following the same reasoning wc find the constraint 



l'5o| 



< 



MX {uJemn — ^I'm'n') 



M . 



(21) 



Consider for example the relative amplitude between the fundamental mode and the first overtone. For Schwarzschild 
black holes and small values of n the typical difference in the imaginary part of the frequency for two consecutive 
overtones (£' = m' = m, n' = n + 1) is 



Ml {ujfr, 



~ 0.2. 



Setting again e = 10~^ the maximum allowed uncertainty on the starting time would be quite small: 5o < 0.05Af 
(this presumably already precludes assuming that the perturbation propagates with speed one, as in flat spacetime). 

Suppose we want to resolve corotating and counterrotating components of the fundamental mode with 1 = 2 (say, 
the components with m ~ ±£). In the case of a spinning black hole background these QNM frequencies arc different, 
but their imaginary parts are actually quite close for most values of the rotation rate jS. i29i| . For example, looking 
at Table m we see that for spin j = 0.5 the difference is \MI {10220 ^^^2-20) I — 0.00212, so that (5o ^ 4.7M. Even 
for a rapidly rotating black hole with j = 0.9 the difference is not as large as between a fundamental mode and its 
overtone: MI {uj 220 - W2-20) — 0.0245, and 5o < 0AM. 




Figure 4: Critical uncertainty in the starting time, as defined by Eq. H21|l . assuming e = 10~^. In the left panel we give the 
critical So for fundamental modes {n = n' = 0) with different angular dependence. For the first mode we assume ^ = m = 2; 
the second mode has I' = 2 and different values of m' = 1, 0, —1, —2 (lines from top to bottom). In the right panel we show 
the critical uncertainty in the relative amplitude of the fundamental mode and first overtone, i.e., n = and n = 1. Here we 
set I = I' = 2, consider all values of m = m' and once again we assume e = 10~^. 



Critical starting-time uncertainties for e = 10~^, general values of the spin and different pairs of modes are plotted 
in Fig. 0] Determining the relative amplitude of a fundamental mode and of the first overtone is generally harder, 
unless we consider corotating modes and near-extremal black holes, as we do in Sec. IIII Fl The spin dependence of i5o 
is quite weak for overtones, but (5o can change by orders of magnitude for modes with different angular dependence 
(t ^ i' ox m ^ m'). For j < 0.5 the time-shift problem, as we defined it here, becomes irrelevant when we want 
to determine the relative amplitude of components with the same / and different m's. The reason is simply that 
modes with different m's have the same QNM frequency in the Schwarzschild limit, so that 5q ^ 00. As a rule of 
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thumb, determining the relative amphtude of angular components with the same I and different m's is harder for large 
rotation. However, as we said before, even for j = 0.9 the critical uncertainty is Sq > 0AM , an order of magnitude 
larger than the typical uncertainty to resolve overtones (which in most cases is ~ 0.05 A/). Most of the qualitative 
features of Fig. 0] are also seen in the experimental problem of resolving different QNMs in the actual detection of a 
ringdown signal (compare e.g. Figs. 3, 4 and 18 of Q). 

In Sec. lIIIEl and Sec. lIIIFl wc will study in more detail the extraction of corotating and counterrotating modes and 
of overtones, respectively. In preparation for this study, in the next Section we outline the general method by which 
we extract quasinormal frequencies from our numerical waveforms. 



C. Extraction of QNM frequencies by an optimal choice of the fitting interval 

Once we have the different multipole components of the numerical solution, we analyze them by applying a fitting 
procedure to each of these components. Since each mode decays exponentially while oscillating with its quasinormal 
frequency, the obvious function to fit the numerical waveform is Eq. (|16|) . where the free parameters are the amplitudes, 
phases and frequencies. As discussed in Sec. IIII Fl only in some cases we have been able to fit for overtones, in the 
sense of getting their expected quasinormal frequencies with reasonable accuracy. However, as described below, the 
residual that we get by truncating the sum at the fundamental mode is already quite small (see also Fig. El • 

In this subsection we are interested in extracting the quasinormal frequencies from our numerical data. To a very 
good approximation the frequencies are independent of to, and we can therefore pick any value for the latter. We still 
need to find a good choice for the time interval [Ti, Tf] over which the ringdown dominates and the fitting procedure 
works best. Since in principle the parameters obtained from the fitting might depend on the choice of this time 
interval, we discuss our procedure in detail. 

Only during the ringdown phase does the waveform have the functional behavior of Eq. (|16() , so the time interval 
[Ti,Tf] should not include the transient regime and the tail phase. For our simulations we found it reasonable to 
pick Tf = 150A'/, since for T > Tf the system typically goes into the tail phase. The choice of T,; is more delicate: 
small values would bring the fitting time window out of the ringdown phase, but large values would make the fitting 
interval small and the resulting fit inaccurate. We decided to take a pragmatic approach: for different values of Ti we 
compute the (relative) residual R{Ti,to) between the fitted function and the numerical data, which we define as 




We then choose the value of Ti that minimizes the residual. In a very well defined sense, this gives an optimal choice 
for Ti. In principle one could use other norms (for example, a sum over squares instead of a sum over absolute values), 
but we checked that this does not affect significantly the results of this paper. Choosing the value of Ti that minimizes 
the residual defined above should not be confused with the minimization procedure done at each Ti to get the fit 
itself. 

Instead of extracting the quasinormal frequencies through a fitting procedure, in principle one could also perform a 
Fourier transform of the solution, as in Ref. |2(il |. However we have found that the fitting procedure provides us with 
far superior accuracy, even in cases with relatively few sampling points. Nonetheless we compared to the results that 
we obtained by Fourier analysis and found consistency between both methods. 

Figure [S] shows the residual as a function of Ti for one of our simulations (the one corresponding to spin j = 0.5 
and i = m = 2 initial data in Table OJ. The residual is independent of the choice of excitation time to, since a change 
in to is just absorbed in the amplitude of the fitting function, leaving the other fitting quantities unaffected. 

Since the black hole's spin is non zero, both m = 2 and m ~ —2 modes are present in the solution. Here we discuss 
only the m = 2 part of the numerical solution. The m = —2 part behaves similarly (in Sec. IIII El we present a detailed 
study of the relative amplitudes of corotating and counterrotating modes). 

From Fig. El we see that i?(Ti,to) has a rather sharp local (and global) minimum. By computing the derivative 
(through finite differences) of the residual with respect to Ti wc find that the minimum is located at Ti — (59.65 ± 
0.025)M. The uncertainty refers to the difference between two consecutive values of Ti, which is in turn given by the 
time step for this simulation: At ~ 0.025M. 

Figure El shows the real and imaginary parts of the frequency extracted from the same simulation as a function of 
Ti. By evaluating them at T^ = (59.65 ± 0.025)M we get lur = 0.585887 ± 1 x 10"^ and ivj = 0.0933851 ± 5 x 10"'^. 

Figure El also reveals that uj changes very little within the interval 50M < Ti < 80M. Since our choice of Ti is 
by no means unique — for example a different definition of the residual would slightly shift Ti — this plateau in the 
frequencies guarantees that the physical quantities we extract are not too sensitive to that uncertainty. This means 
that the errors in our numerically extracted QNM frequencies due to the choice of Ti are quite small. 
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Figure 5: Residual in the fit, as defined in Eq. H22|l . as a function of the initial time for the fitting Ti. Looking at the minimum 
of the residual we can determine Ti with high precision. This plot corresponds to a simulation with spin j — 0.5, I = m — 2 
initial data with A — 0, B = 1 and a radial resolution Ar — M/10. 
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Figure 6; The left and right panels show the real and imaginary parts of the quasinormal frequencies extracted from the 
simulations in Fig. |S] From the optimal starting time determined by minimizing the residual, Ti = (59.65 ± 0.025)M (see 
previous figure), we find ujr = 0.585887 ± 1 x 10"" and loi = 0.0933851 ± 5 x 10"^ 



We are now ready to examine the quasinormal fre quen cies obtained from our numerical data in the way just 
described. Table shows the frequencies computed in [13] using Leaver's continued fraction method for perturbed 
Kerr black holes with spin j = 0, 0.5, and 0.9 (here labeled perturb.). Along with these frequencies we list values 
extracted from our time domain evolutions (labeled numerical) and the relative differences between the two. The 
numerical values were obtained by evolving different initial data sets with A = and (£ = 2, m = 0, ±2) in Eq. H10|l . 
and fitting for the multipoles present in the initial data (we discuss the additional multipolcs generated by rotational 
mode mixing below). For j = the frequencies do not depend on m, therefore we only show results for to = 0. Even 
with a relatively modest resolution, the differences on quasinormal frequencies from our three-dimensional simulations 
in Table ^ arc between one and two orders of magnitude smaller than the ones reported in previous two-dimensional, 
axisymmctric simulations of gravitational perturbations . 



D. Rotational mode mixing 

In Sec. IIIDl we described our initial data family sets, which were expanded in spherical harmonics. Since the Kerr 
background is not spherically symmetric we should not expand the perturbation in terms of spherical harmonics, 
but (more rigorously) in terms of the spin- weighted spheroidal harmonics s^imio-^)-, where s is the spin weight of 
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the perturbing field, a = jM is the black hole spin parameter, and uj is the frequency in a Fourier expansion of the 
perturbation (a quantitative discussion of spin-weighted spheroidal harmonics and more references can be found in 
|3flj| ) . However, as first shown by Press and Teukolsky the sSem^s may be expanded as a power series in auj: 

sSim = sYim + {auj) ^ ctim sYtm + 0(aw)^ . (23) 

Here sYim denotes a spin-weighted spherical harmonic of spin-weight s. In this paper we focus on scalar perturbations 
(s = 0), in which case the spin- weighted spherical harmonics reduce to ordinary spherical harmonics. The coefficients 
ci'im are related to the more familiar Clebsch-Gordan coefficients [s^,!^- As a result of (|^ . and because of the 
orthogonality of the (spin-weighted) spherical harmonics, inner products of different spheroidal harmonics will be 
given by inner products of spherical harmonics with higher-order corrections in auj. At least for small aoj, we may 
expect these contributions to be small. In fact, the corrections turn out to be small even for moderately large values 
of au (see for an explicit calculation of the inner products at the QNM frequencies). Nevertheless, using spherical 
harmonics instead of spheroidal harmonics can induce a small amount of mode-mixing in the initial data. 

For a spherically symmetric background spacetime, initial data with different values of t evolve separately and the 
angular structure of each mode is preserved during evolution. On the other hand, for a Kerr background with nonzero 
spin, modes with different values of £ do couple and furthermore, modes that are not present in the initial data can be 
excited during evolution. This may make it necessary to increase the angular resolution compared to the non-rotating 
case to resolve the higher i modes generated during evolution. However, the decay rate of these modes increases with 
I, so even when modes with higher values of (. are generated during evolution, they do not dominate. Therefore, we 
found that if we accurately resolve the angular structure initially, the same is in general true for the whole evolution. 

Figure [3 illustrates rotational mode coupling for non-zero spin backgrounds (see also 1^3 arid |33l | for numerical 
studies of mode-mode coupling). Since modes with same m but different £ can couple to each other, we show the 
extracted (£ = 4, m = 2, n = 0) waveform (for three simulations with different spin parameters) excited by initial 
data whose angular dependence is given by an £ = to = 2 spherical harmonic. As expected, the rotationally-induccd 
excitation of the (£ = 4, to = 2) mode typically increases with spin. Some additional mode mixing is an artifact of 
the symmetry of our computational grid. This "spurious" mode mixing is present also for j ~ 0, but it converges to 
zero as we increase the angular resolution. All other modes we extract, up to £ = 4 and all allowed values of to, are 
within roundoff error throughout the simulations. 
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Figure 7: The left panel shows the extracted {£ = 4,m — 2,n = 0) waveform for three simulations with different spin parameters 
as seen by an observer at r = 5M. The initial data are a pure {£ — m = 2) mode and are set up according to Eq. llOt with 
A = 0, B = 1 and ro = 20M. For zero spin the different multipole components of the solution should evolve independently 
and no modes besides the one in the initial data should be excited, while for non-zero spin modes with different £ but same m 
do couple [2^ . In the Schwarzschild case the (£ = 4, m = 2, n = 0) waveform differs from zero due to our grid structure and 
discretization errors, but it converges to zero with increasing resolution. This is illustrated by the right panel, which shows 
the extracted (£ = 4, m = 2, n = 0) amplitude for j = and j — 0.9 from runs with two resolutions (20 x 20 x 1000 and 
30 X 30 X 1500 points per patch and outer boundaries at lOOM). Only for j = 0.0 the mode converges to zero. 



Since we only extract QNMs up to £ = 4 we need to test whether there is a relevant contribution from higher modes 
that we do not extract explicitly. In the absence of £ > 4 modes, summing up all extracted modes up to £ = 4 we 
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should recover the full field, according to Eq. (|12|l . The result of this test for a spinning black hole with j = 0.9 is 
shown in Fig.|Sl at the level of accuracy needed in the present work, extracting modes with ^ < 4 is sufficient. 
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Figure 8: Results from a run with initial data parameters £ = m = 2 and spin j = 0.9. The left panel shows the square of the 
amplitude of all modes up to ^ = 4 which are not within the roundoff error. The right panel shows the sum over the square of 
those modes compared to the integral over a sphere of the full field squared: As expressed in eauation ll2l the two curves lie on 
top of each other, and there is no relevant contribution from higher modes. 



E. Relative amplitude of corotating and counterrotating modes 

We know that when the solution is in the quasinornial ringing regime, it will behave according to Eq. Hlt)|) . In the 
previous subsection we have verified through our simulations the values predicted in Ref. [l6| for the frequencies. We 
now also want to verify the amplitudes of each mode, as predicted in that same reference. 

Assume that the observer and the initial data are located far away from the black hole (these assumptions underlie 
the "asymptotic approximation" adopted in From Eq. (4.15) of when B = the response of the black 

hole in the ringdown phase should be well approximated by a QNM decomposition of the form 

<^im{r,t) « -^^/^(77^ |x](iAu;f™„)i?£™„e-'^'"?""/4g-i^,„„(t-ro-r.)| ^ (24) 
In our simulations we set A = 0, in which case it can easily be shown that the previous expression becomes 

^imir, t) « ^^V^aBU I i?^™„e-^'"-""/'^e-'"-"" I , (25) 

^ [ n=0 J 

With respect to we added an extra factor ro/r. This is because Eq. (4.15) in refers to the Sasaki-Nakamura 
function X^^_^ (r, t) , which is related to the Teukolsky function $f „ ('', t) that we are using in our evolutions by the 
relation xj^^(r,t) = (r^ + a^) $^„(r, <) (see the discussion in Appendix C of |l6l|l. We are interested in large 

values of r, for which the asymptotic approximation holds and x£"^'(r, t) ~ r^im(r,t). The transformation between 
the Teukolsk y a nd Sasaki-Nakamura functions must also be taken into account when comparing the initial data in 
Eq. (4.14) of (Tg with our initial data, Eq. (|10f) . Assuming a vq and r 3> 1 this comparison yields the normalization 
factor rp in the equations above. 

The scalar QNM frequencies ujimn and the scalar excitation factors Bimn are listed in Table I and Table III of 
|l6j . respectively. In that reference and in Eq. (|24|l Boyer-Lindquist coordinates are used; since in our simulations we 
use Kerr-Schild coordinates we need to transform Eq. H24|l appropriately. Since $ is a scalar, the transformation is 
straightforward. The transformation of the initial data is more subtle, since the slices are different. One would expect 
that whenever the asymptotic approximation is valid the difference between the slices should not be too important. 
The results discussed below and explicit comparisons between evolutions using both coordinate systems in the non- 
spinning case jssf confirm this expectation. Details on how we transform the initial data and the field itself are given 
in Appendix fXl 



14 



To check the accuracy of Eq. (|25() . in the rest of this section we analyze evolutions of different initial data sets, 
all of them consisting of a combination of = 2, m = 2) and (£ = 2,to = —2) modes with A = and B = 1. We 
numerically explore the dependence of the amplitudes of the counter- and co-rotating fundamental modes (in the 
next subsection we will study overtones) on the width a of the initial data [cf. Eq. H1U|) ]. In order to assess more 
quantitatively the effect of the time-shift problem (see Sec. IIII B|) wc first compare the value of the width maximizing 
these amplitudes. Given that all the initial data sets that we consider are centered at the same radius, we can 
make the reasonable assumption that locally (that is. around the width for which the amplitudes are maximal) to is 
approximately the same for each set. If to were exactly the same, the value of to used would not change the width at 
which the maximum amplitude is located, since changes in to would only involve a global rescaling of all amplitudes, 
as discussed in Sec. IIII Bl Therefore the hope is that within the setting described for our simulations the width for 
which the amplitudes are maximal does not depend too sensitively on to. 

The numerical results shown here were obtained with the same number of points in the angular direction as above. 
We used half the resolution in the radial direction (that is, Ar = M/5) for a rough scan of a large a range, and again 
the original resolution around the maxima of the amplitudes. We chose initial data with varying widths a, rg = 20M 
(as in the simulations above) and an observer at r = 40Af, for which the asymptotic approximation holds reasonably 
well [3^ . We picked to = Hnitiai data + f (that is, to = 60M in the cases considered), which is approximately the time 
the initial data pulse needs to propagate to the black hole and back to the observer. 
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Figure 9: Numerically obtained excitation amplitudes of the £ — \rn\ = 2 fundamental modes assuming an observer location 
'"obs = 40M and a ringdown starting time to = 60M. The left panel refers to a black hole with spin j = 0.5. According to 
predictions from perturbation theory in the asymptotic approximation [cf. Eq. 12411 and the following discussion] the maximum 
for m — 2 should be located at (J220 ~ 2.445, while the value that we obtain from our simulations is 0-220 ~ 2.55 ± 0.05 
(the uncertainty describing the difference between consecutive values of a used in our simulations: Act = 0.05). Similarly, for 
m = —2 the width at the maximum should be 0-2-20 = 3.434, while we obtain 02-20 = 3.875 ± 0.075. The right panel, in 
turn, refers to a black hole with spin j = 0.9. In this case the theoretical (numerical) maxima are located at 0220 = 1.816 
(0220 = 1.85 ± 0.05) and 02-20 = 3.758 (02-20 = 3.85 ± 0.05), respectively. The inset in the left panel is a zoom around the 
maximum for j — 0.5 and m = 2. As discussed in the text, an uncertainty in the excitation time of 0.09M would already 
explain the difference between the predicted location of the maxima and our numerical results. 



Figure|51and (more quantitatively) Tables UTI and IIIII show the excitation amplitudes as functions of the width of the 
Gaussian a from our numerical simulations. At first our results could be interpreted as an approximate verification 
of the predictions of |l6l |. However, if one takes into account the limitations imposed by the time-shift problem, 
the agreement can in fact be considered excellent. For example, take the j = 0.5, m ~ —2 case, which is the one 
where the difference between the theoretical and numerical values is largest. The theoretical maximum is located at 
a = 3.434A'/, while the numerical value is cr = (3.875 ± 0.075)M (the uncertainty indicating the difference between 
consecutive values of a). The relative numerical amplitude between cr = 3.45M and cr = 3.85A/ — 3.9M from our 
simulations is « 1.008 (see Table lll|l . If the values of to for these two widths differ by « 0.09A'/, the amplitude 
corresponding to cr = 3.45Af would actually be larger than the one of cr 3.85Af — 3.9Af and would therefore shift 
the maximum to the predicted value of 3.45Af . Recalling that we used to = 60Af, a very modest uncertainty in the 
relative ringdown starting time (« 0.4%) would shift the maximum to the theoretical value. We also assumed the same 
excitation time to for all the initial data sets when fitting our numerical data. Whenever such assumption is a good 
approximation, the precise value of to should not affect the location of the width for which the excitation amplitudes 
is maximal. In particular, the approximation should be good if the initial data pulses are relatively narrow. However, 
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Table II: Excitation amplitudes for j = 0.5, £ = 2 and n — for initial perturbations of variable Gaussian width a, as displayed 
in Fig. |5] The observer location in these runs is r = 40M. Highlighted are the maxima in the amplitudes of the different 
m-modes. Also shown are the relative amplitudes of the two modes, and the relative differences between the values predicted by 
perturbation theory and the ones extracted from our numerical simulations. The amplitudes are given for the wave expressed 
in Boyer Lindquist coordinates (see appendix EI for details) and are multiplied by a factor of r/ro to get them in an observer 
independent form. 
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as a increases, the possibility of the excitation time to shifting around has to be taken into account, because the 
interaction time of the pulse with the black hole becomes longer and the interaction sets in well before the center 
of the pulse reaches the black hole. Taking all this into account, the agreement between numerical and perturbative 
results for the location of the maxima can be considered excellent. The situation for the amplitudes themselves is 
different, as discussed next. 

Tables m and IIIII show the predicted and extracted absolute and relative amplitudes for the co- and counter- 
rotating modes Am=2, ■Am=~-2- As expected, the prediction from perturbation theory works better for sharp pulses. 
The differences between the predicted and absolute values are of order a few percent for sharp pulses and grow with 
cr. For (7 = 4 the difference is as large as ~ 20% and ~ 60% for j = 0.5 and j = 0.9, respectively (the actual 
amplitudes being larger than the predicted ones) . These large differences in the relative amplitudes are mostly due to 
the amplitude of the corotating mode, the predicted and extracted amplitudes for the counterrotating one agree quite 
well. The fact that the location of the maxima, as discussed above, agrees very well despite the large differences in 
the amplitudes for larger tr can be easily explained: the location of the maxima for the corotating mode takes place 
at tr « 1.85Af , which corresponds to a pulse which is sharp enough for perturbation theory to give a good prediction, 
while the maximum for the counterrotating mode is at a larger value of a but, as we have discussed, the agreement 
between predicted and measured amplitudes is quite good for that mode. 

Could this large difference in amplitudes be explained by the time-shift problem, as discussed in Sec. IIII Bf Using 
Eq. (|20|) and assuming that t^ is roughly the same for both modes we find that an uncertainty in the excitation time 
as large as 5q = ±5A/ would hiiply an uncertainty on the relative amplitudes of about ±1.1% for j ~ 0.5, and ±13% 
for j = 0.9. Therefore the uncertainty 5^ does not seem to account for the differences that we find with respect to 
the predicted amplitudes. One possibility is that the excitation time to is different for the two modes in a pair; but, 
if so, it is not clear then why our naive choice of to is very good for the counterrotating mode and quite bad for the 
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Table III: Same as TableEllfor j = 0.9. 





numerical results 


perturbation theory 


relative difference 


a 


.4220 


.42-20 


.42-2o/.4220 


.4220 


.42-20 


.42-2o/.4220 


.42-20 /.4220 


1.60 


0.1594 


0.3156 


1.98 


0.1768 


3683 


2.08 


0.05 


1.70 


0.1615 


0.3319 


2.06 


0.1766 


0.3857 


2.18 


0.06 


1.75 


0.1621 


3399 


2.10 


0.1755 


3990 


2.27 


0.08 


1.80 


0.1625 


0.3476 


2.14 


0.1752 


0.4022 


2.30 


0.07 


1.85 


0.1626 


0.3553 


2.18 


0.1740 


0.4101 


2.36 


0.07 


1.90 


0.1625 


3629 


2.23 


0.1725 


0.4177 


2.42 


0.08 


2.00 


0.1617 


0.3775 


2.33 


0.1725 


0.4323 


2.51 


0.07 


3.60 


0.0800 


0.5173 


6.47 


0.0566 


0.5259 


9.30 


0.30 


3.70 


0.0743 


0.5192 


6.99 


0.0507 


0.5235 


10.32 


0.32 


3.75 


0.0714 


0.5204 


7.29 


0.0479 


0.5220 


10.89 


0.33 


3.80 


0.0688 


0.5204 


7.56 


0.0452 


0.5203 


11.51 


0.34 


3.85 


0.0661 


0.5212 


7.89 


0.0427 


0.5184 


12.15 


0.35 


3.90 


0.0636 


0.5208 


8.19 


0.0402 


0.5162 


12.85 


0.36 


4.00 


0.0590 


0.5194 


8.81 


0.0355 


0.5114 


14.41 


0.39 


4.10 


0.0590 


0.5184 


8.79 


0.0312 


0.5059 


16.20 


0.46 


4.20 


0.0505 


0.5144 


10.19 


0.0274 


0.4997 


18.25 


0.44 


4.30 


0.0469 


0.5098 


10.88 


0.0239 


0.4929 


20.64 


0.47 


4.40 


0.0433 


0.5047 


11.65 


0.0207 


0.4854 


23.41 


0.50 


4.50 


0.0433 


0.4991 


11.52 


0.0179 


0.4774 


26.63 


0.57 



corotating one. It is actually not clear why such a large disagreement happens only for the corotating mode, and not 
for the counterrotating one. The possibilities that the initial data and/or the observer are not far enough away for 
the asymptotic approximation to be valid, or that the d isag reement is due to a lack of resolution, seem to be ruled 
out by one-dimensional studies in the non-spinning case [33 ■ Summarizing, even though the exact mechanism is not 
clear, all this suggests that the predicted amplitudes for the corotating mode in the asymptotic approximation are 
simply valid only for very sharp pulses, as the black hole spin increases. 

To conclude, we want to discuss one aspect of our simulations, as shown in Figure El We see a rather large 
discrepancy between the amplitudes resulting from runs with resolution A?' = il//5 and Ar = M/10, especially for 
j = 0.9. That is a direct effect of decreasing accuracy in X(u;£m„) when going to high spins (see Sec. IIII C|) and the 
need for more resolution in those cases. The location of the maximum, however, is always consistent (that is, within 
the differences in a used in the different initial data sets) between runs of different resolution. That is not surprising 
since the measured w^mn at a fixed resolution is roughly the same for all values of tr, and the value of a that maximizes 
A only depends on the value of cuemn- 

F. Overtones and rapidly spinning black holes 

As discussed in the introduction, a single complex quasinormal frequency contains enough information to determine 
the two parameters of a Kerr black hole (namely, its mass M and spin j). If one is able to detect a second mode 
from the same source, one can use this extra information for a consistency check that would increase the confidence 
in the interpretation of the measured data as signals from a perturbed black hole. An important question that might 
be answered by numerical relativity is whether more than one mode will be detectable by Earth- and space-based 
gravitational wave detectors. In Sec. IIII El we considered the relative amplitude of corotating and counterrotating 
modes; here we use our simulations to determine the relative excitation of overtones with the same angular dependence 
and m > 0. According to perturbation theory, in this case the damping time of the first overtone becomes comparable 
to the damping time of the fundamental mode for large spins (see Fig.0}. In addition, the excitation factor of higher 
overtones is usually larger than the excitation factor of the fundamental mode for large j |l6j | . This means that higher 
overtones are more likely to be detectable for fast spinning black holes. A detailed study of this topic is beyond the 
scope of this paper, but here we briefly discuss how we can extract information about overtones from our data and 
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determine which modes contribute most significantly to the waveform. 

We perform simulations for different spins {j = 0, 0.5, 0.9, and 0.98). The initial data and numerical procedure are 
the same as in Sec. Ill D1 and III El with one exception: for spins j > 0.9 wc found it necessary to increase the angular 
resolution. The simulations presented in this section used a resolution of 31 x 31 grid points on each block in the 
angular directions. This is not surprising, since for fast rotation we expect more dynamics in the angular directions. 

The extraction of modes is done in principle according to Sec. lIII C\ Extracting information about all modes present 
in the data can turn into a subtle problem, especially when the contributions of some modes is weak. One option 
is to first fit for the strongest mode present in the data, subtract the fit, fit for the next dominant mode and so on, 
repeating the procedure as long as an oscillatory exponential decay is seen in the data. However, when there are 
several modes with similar contributions we can just fit for all of them at the same time. This is exactly what we did 
for fundamental modes with different m in the previous subsection. When a single mode dominates the waveform the 
first strategy not only seems to be more meaningful, but also turns out to work better in practice. The results of this 
section were computed by a hybrid of these two methods, depending on the contribution of each mode (something 
that one can find out by, for example, looking at the dominant frequencies of the signal to fit). 

Table llVl shows the quasinormal frequencies of the overtones that we get from our simulations, using (A ~ 0, B = 1), 
a = A/, ro = 20M and an observer at r = 60A/. We find that the overtones for the m = —2 mode do not contribute 
enough to the waveforms to extract them from our data with decent accuracy, especially for high spins. The reason 
for this is that the imaginary part of their frequency is generally smaller than the one for the corresponding m ~ 2 
mode, which makes them decay faster. The decay of the m = 2 mode, on the other hand, slows down considerably 
when increasing the spin. We numerically find that the excitation amplitude (at fixed to) increases with increasing 
spin. Those two effects combined make the extraction of overtones easier and more accurate in the high spin cases. 
Quite remarkably, for runs with spin j = 0.9 and above we can extract the quasinormal frequency for n = 2 with 
reasonable accuracy (see Table llV)l . 

Table IVl compares the amplitudes of the three most dominant / = 2 modes, (to — 2,7i ~ 0), (to = — 2,n = 0) and 
(to = 2,n = 1), with the predicted asymptotic amplitudes of Eq. H24(l . Except for the j = 0.98 case, the difference 
between the predicted and extracted values for the relative amplitudes between a given mode and the fundamental 
^ = 2 = m one is of the order of a few percent for the fundamental mode and one order of magnitude larger for the 
first overtone. 

Table IV; Comparison of quasinormal frequencies for the first overtones (n = 1, 2) of an = 2, m = 2 mode, for black holes with 
varying spin, as predicted by perturbation theory and as extracted from our numerical simulations, along with their relative 
differences. This table is complementary to Table where we show the frequencies associated to the fundamental modes. The 
extraction of overtones becomes easier for rapidly rotating black holes, as explained in the text, allowing us to extract the 
frequencies of two overtones for high spins. 



j 


n 


'-^-'perturb 


^numerical 


rel. difference (Re, Im) 


0.0 


1 


0.46385 - 0.29560i 


0.45651 - 0.288591 


1.58 X 10"^ 2.37 X 10"^ 


0.5 


1 


0.57344 - 0.283341 


0.54718 - 0.317221 


4.58 X 10"^ 1.20 X 10"^ 


0.9 
0.9 


1 
2 


0.77768 - 0.20801i 
0.77043 - 0.34720i 


0.73737 - 0.195581 
0.52473 - 0.353191 


5.18 X 10"^ 5.98 X 10"^ 

3.19 X 10"\ 1.73 X 10"^ 


0.98 
0.98 


1 
2 


0.89622 - 0.12214i 
0.89358 - 0.202441 


0.93152 - 0.124061 
0.88668 - 0.258501 


3.94 X 10"^ 1.57 X 10"^ 
7.72 X 10-^ 2.77 x lO'^ 



IV. CONCLUSIONS 

The chances of a multi-mode detection by either Earth- or space-based gravitational wave detectors will depend on 
the relative amplitude of those modes. Knowing in advance which modes should be excited under a realistic binary 
merger would reduce the dimensionality of the template bank needed to perform matched filtering on ringdown 
waveforms. An answer that numerical relativity might provide is precisely which modes are likely to be dominant. 
This involves predicting the relative amplitudes of different pairs of modes under a variety of scenarios. In this paper 
we have taken a first step towards understanding the issues involved in such a prediction. 

We first presented a systematic way of extracting QNMs from a given signal. Our procedure has a number of built 
in self-consistency checks, to make sure that when we keep adding modes to our fit we are fitting a true signal and 
not numerical noise. One of these self-consistency checks is to make sure that we extract the correct quasinormal 
frequency of each mode within a certain accuracy. If the data being analyzed comes from a numerical simulation. 
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Table V: Absolute and relative amplitudes as a function of the black hole spin and angular dependence of the perturbations, 
as predicted by perturbative calculations and as extracted from our numerical evolutions. The amplitudes are given for the 
wave expressed in Boyer Lindquist coordinates (see appendix^for details) and are multiplied by a factor of r/ro to get them 
in an observer independent form. The last column presents the relative difference between perturbative and numerical results 
for relative amplitudes. In the corotating case we also extract the amplitude of the first overtone. The differences in the 
relative amplitudes are considerably smaller when we look at corotating and counterrotating modes, compared to the case of 
fundamental mode and first overtone with the same angular dependence. This can be explained by the relative magnitude of 
their damping frequencies, as discussed in Section fill Bl fsee also Table HVl . This difference becomes less pronounced at very 
large spins, as expected from the analysis of Section Till Bl 





mode 


numerical result 


perturbation theory 


relative difference 
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I m n 












0.00 


2 2 


0.211 


1.00 


0.221 


1.00 


0.00 


0.00 


2 2 1 


0.316 


1.50 


0.504 


2.28 


0.342 


0.50 


2 2 


0.201 


1.00 


0.213 


1.00 


0.00 


0.50 


2-2 


0.208 


1.03 


0.228 


1.07 


0.037 


0.50 


2 2 1 


0.525 


2.61 


0.768 


3.61 


0.277 


0.90 


2 2 


0.137 


1.00 


0.148 


1.00 


0.00 


0.90 


2-2 


0.211 


1.54 


0.246 


1.66 


0.072 


0.90 


2 2 1 


0.533 


3.89 


0.98 


6.62 


0.412 


0.98 


2 2 


0.0833 


1.00 


0.068 


1.00 


0.00 


0.98 


2-2 


0.263 


3.16 


0.257 


3.78 


0.164 


0.98 


2 2 1 


0.634 


7.61 


0.416 


6.12 


0.243 



consistent frequencies can be used to monitor the accuracy of the code. If the data is experimental, consistency of 
the frequencies allows for a test of the no hair conjecture. In more detail: during our fitting procedure we first fit for 
the dominant mode(s), look at the residual (defined as the difference between our original signal and the fit), make 
sure that it has a consistent quasinormal ringing behavior and only then fit for the next set of modes, repeating the 
procedure as long as it makes sense to do so. By following this procedure we have hardly been able to go beyond 
the first few dominant modes, and this was only possible in very special cases. We expect this to happen with most 
numerical simulations. 

We addressed in some detail the so-called time-shift problem. In essence, this is the fact that the quasinormal 
amplitudes depend exponentially on the quasinormal ringing excitation time, which is not defined unambiguously 
(not even in the continuum). Furthermore, examining actual values of quasinormal frequencies we have seen that 
this exponential dependence is an important factor to take into account in practice. To (partially) get rid of this 
exponential dependence wc propose to look at relative amplitudes: choosing pairs of modes whose damping frequency is 
as close as possible, we can partially cancel each others exponential dependence. We analyzed in detail the exponential 
dependence of different pairs of modes as a function of the black hole spin. In particular, we found that the time-shift 
problem becomes more important as one increases the spin. For modes with the same value of £, for example, the 
problem is not very relevant for spins j < 0.5. On the other hand, an accurate extraction of the relative amplitude 
between the fundamental mode and the first overtone only seems feasible for very high spins and m > 0. 

Keeping this in mind, we first extracted the fundamental quasinormal frequencies for different values of spin, ranging 
from J = to a rapidly rotating black hole with j = 0.98. Even using modest resolutions our frequencies agree with 
those obtained from perturbation theory within one part in 10^ to one part in 10^, depending on the black hole spin, 
location of the observer and angular dependence. To our knowledge this is the first time that quasinormal frequencies 
for scalar perturbations of Kerr, as predicted by perturbation theory, have been verified by numerical evolutions of 
the field equations. 

Next wc analyzed in detail the relative amplitude of corotating and counterrotating fundamental modes, as a function 
of the width of the initial perturbation and the black hole spin, being able to quantify (within the limitations imposed 
by the time-shift problem) under what conditions the asymptotic approximation of Rcf. |l6l | is valid. In particular, 
we were able to verify the widths of the initial perturbation corresponding to the maximal QNM excitation. Finally, 
we studied the excitation of overtones. Wc found that, according to expectations from perturbation theory |l6l |. they 
get significantly excited for corotating modes and very high spins. In this particular case we were able to extract the 
complex QNM frequency for the fundamental mode and the first two overtones, with a difference with respect to the 
predicted values by perturbation theory of the order of a tenth of a percent to ten percent, depending on the mode 
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and the black hole spin. We expect the techniques and results of this paper to be general enough to be useful for 
future work on ringdown waveforms. 
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Appendix A: CHANGE OF COORDINATES AND INITIAL DATA 



To compare the numerical results with predictions from perturbation theory we must switch from the usual Boyer- 
Linquist coordinates (as used, for example, in to the Kerr-Schild coordinates used in our code. We denote 

by (r^, t) the Boyer-Lindquist radial tortoise coordinate and time, and by (f, t) the Kerr-Schild coordinates. The 
transformation we need is given by 

t{t,f) = t-n{f)+i, (Al) 

n(f) = f + n{f), (A2) 

with the definitions 

r+ = M + - a2 , (A4) 

r_ = M - \Jm^ - a2 . (A5) 

(A6) 

The reference time t can in principle be freely chosen and is used to define where t(f , i) crosses zero. We fix it by the 
condition that in both coordinate systems the initial pulse is at the same physical distance from the black hole, i.e. 
i(i = 0,f = fo) = 0: 

f=r!(fo). (A7) 
The location of the initial pulse in these coordinates becomes 

ro = fo + n{fo). (A8) 
For consistency, the value of a has to be adjusted to tortoise coordinates. As a rough approximation we set 

(T ^]^[r*{f + a)-r^{f -a)]. (A9) 

From eauations ll6l and l25l we read of the exponential decay of each modes amplitude in Boyer-Lindquist coordinates 
and then substitute them by the Kerr-Schild coordinates. We choose = + ?'o- 

This equation relates the amplitudes as seen in Boyer-Lindquist coordinates Aemn with the ones found in the simu- 
lations that were done in Kerr-Schild coordinates Aimn- 
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